% fig2c.m ——— 复刻论文 Fig.2(c)：平均哈密顿能量曲线 ———
clear; clc; close all;

%% —— 全局参数 —— 
beta   = 1;      % 论文中 β
xi     = 0.8;    % 论文中 ξ
mu     = 1;      % 论文中 μ
A      = 3.4;    % 外激励幅值 A
omega  = 1;      % 外激励角频率 ω

dt     = 0.01;   % 时间步长
T0     = 400;   % 丢弃瞬态时间
Tsim   = 600;   % 统计时间
N0     = floor(T0/dt);
Nsim   = floor(Tsim/dt);

S0     = [0.2; 0.1; 0.2];   % 初始状态 [x; y; z]

%% —— 扫描 α 范围 —— 
alpha_min = 1.0;
alpha_max = 4.0;
Nalpha    = 200;
alpha_list = linspace(alpha_min, alpha_max, Nalpha);

H_avg = zeros(1, Nalpha);

%% —— 并行循环计算平均能量 —— 
parfor ia = 1 : Nalpha
    alpha = alpha_list(ia);
    S = S0;
    Hsum = 0;
    
    % 仿真并累加能量
    for k = 1 : (N0 + Nsim)
        t = (k-1) * dt;
        mu_s = A * cos(omega * t);
        
        % RK4 一步积分主系统
        k1 = f(S, alpha, beta, xi, mu, mu_s);
        k2 = f(S + 0.5*dt*k1, alpha, beta, xi, mu, mu_s);
        k3 = f(S + 0.5*dt*k2, alpha, beta, xi, mu, mu_s);
        k4 = f(S +     dt*k3, alpha, beta, xi, mu, mu_s);
        S  = S + dt*(k1 + 2*k2 + 2*k3 + k4)/6;
        
        % 丢弃瞬态后累加能量
        if k > N0
            xk = S(1);
            yk = S(2);
            Hsum = Hsum + ( xk^2/(2*alpha) + yk^2/2 );
        end
    end
    
    H_avg(ia) = Hsum / Nsim;
end

%% —— 绘图 —— 
figure;
plot(alpha_list, H_avg, 'r-', 'LineWidth',1.5);   % 红色实线
xlabel('\alpha','FontSize',12);
ylabel('\langle H \rangle','FontSize',12);
title('Fig.2(c) 平均哈密顿能量 \langleH\rangle','FontSize',14);
xlim([alpha_min alpha_max]);
grid on;
